/*******************************************************************************
SSC0143 - Programacao Concorrente
Prof. Julio Cesar Estrela

T1 - Calculo de PI com 10 000 000 de casas decimais

Grupo:

Daniel de Oliveira Pereira 6427171
Guilherme Takasawa Yagui 6426698
Vinicius Marques Stocco 6427233

Março de 2013

******************************************************************************** 
Algoritmo de Gauss-Lagendre - Versao Concorrente 
*******************************************************************************/
#include <stdio.h>
#include <math.h>
#include <time.h>
#include <gmp.h>
#include <pthread.h>
#include <semaphore.h>
// numero de casas decimais corretas de PI
#define PREC_DGT 10000000
// numero de bits minimo para armazenar PI com o numero de casas decimais corretas
#define PREC_BIT (int) ceil((log(10) / log(2)) * PREC_DGT)
// numero de iteracoes necessario para calcular PI com o numero de casas decimais corretas desejado
#define ITR (int) ceil(log(PREC_DGT) / log(2))
mpf_t an, an1, an2, bn, bn1, bn2, tn, pn, pi, aux;
/*
calcula a(n+1)
a(n+1) = (a(n) + b(n))/2
*/
void *A(void *args)
{ 
	mpf_add(an2, an1, bn1);	
	mpf_div_ui(an2, an2, 2);
	pthread_exit(NULL);
}
/*
calcula b(n+1)
b(n+1) = sqrt(a(n)*b(n))
*/
void *B(void *args)
{ 
	mpf_mul(bn2, an1, bn1);
	mpf_sqrt(bn2, bn2);
	pthread_exit(NULL);
}
/*
calcula t(n+1)
t(n+1) = t(n) - p(n)*(a(n) - a(n+1))^2
p(n+1) = 2*p(n)
*/
void *T(void *args)
{
	mpf_sub(aux, an, an1);
	mpf_pow_ui(aux, aux, 2);
	mpf_mul(aux, pn, aux);
	mpf_sub(tn, tn, aux);
	mpf_mul_ui(pn, pn, 2);
	pthread_exit(NULL);
}
int main(char argc, char *argv[])
{
	clock_t t; // tempo
	int i, f;
	pthread_t t_a, t_b, t_t, t_p; //threads
	// define a precisão desejada para cada variavel
	mpf_set_default_prec(PREC_BIT);

	// marca o clock inicial da execucao do algoritmo
	t = clock();

	// aloca memoria para as variaveis utilizadas para calcular PI
	mpf_init(an);
	mpf_init(an1);
	mpf_init(an2);
	mpf_init(bn);
	mpf_init(bn1);
	mpf_init(bn2);
	mpf_init(tn);
	mpf_init(pn);
	mpf_init(aux);
	mpf_init(pi);
	// inicializa a(0) = 1
	mpf_set_d(an2, 1.0);
	mpf_set_d(an, 0.0);
	mpf_set_d(an1, 0.0);
	// inicializa b(0) = 1/sqrt(2)
	mpf_set_d(bn2, 1.0);
	mpf_set_d(aux, 2.0);
	mpf_sqrt(aux, aux);
	mpf_div(bn2, bn2, aux);
	mpf_set_d(bn, 0.0);
	mpf_set_d(bn1, 0.0);
	// inicializa t(0) = 1/4
	mpf_set_d(tn, 1.0);
	mpf_div_ui(tn, tn, 4);
	// inicializa p(0) = 1
	mpf_set_d(pn, 1.0);
	mpf_set(an, an1);
	mpf_set(bn, bn1);
	mpf_set(an1, an2);
	mpf_set(bn1, bn2);
	// cria as threads utilizadas para calcular a(n+1) e b(n+1)
	pthread_create(&t_a, NULL, &A, NULL);
	pthread_create(&t_b, NULL, &B, NULL);
	pthread_join(t_a, NULL);
	pthread_join(t_b, NULL);
	for (i = 0; i < ITR-1; i++)
	{
		// calcula pi = ((a(n) + b(n))^2)/(4*t(n))
		mpf_add(pi, an1, bn1);
		mpf_pow_ui(pi, pi, 2);
		mpf_div_ui(pi, pi, 4);
		mpf_div(pi, pi, tn);
		//gmp_printf("Iteracao %d, pi = %.*Ff \n", i, PREC_DGT, pi);
		//gmp_printf("Iteracao %d, pi = %.*Ff \n", i, 100, pi);
		/*
		a(n-1) = a(n)
		b(n-1) = b(n)
		a(n) = a(n+1)
		b(n) = b(n+1)
		*/
		mpf_set(an, an1);
		mpf_set(bn, bn1);
		mpf_set(an1, an2);
		mpf_set(bn1, bn2);
		// calcula a(n+1) e b(n+1)
		pthread_create(&t_a, NULL, &A, NULL);
		pthread_create(&t_b, NULL, &B, NULL);
		pthread_join(t_a, NULL);
		pthread_join(t_b, NULL);
		// calcula t(n+1)
		pthread_create(&t_t, NULL, &T, NULL);
		pthread_join(t_t, NULL);
	}
	// calcula pi = ((a(n) + b(n))^2)/(4*t(n))
	mpf_add(pi, an1, bn1);
	mpf_pow_ui(pi, pi, 2);
	mpf_div_ui(pi, pi, 4);
	mpf_div(pi, pi, tn);
	//gmp_printf("Iteracao %d, pi = %.*Ff \n", i+1, PREC_DGT, pi);
	//printf("\nNumero de iteracoes : %d \n", ITR);

	// libera a memoria alocada para o calculo de PI
	mpf_clear(an);
	mpf_clear(an1);
	mpf_clear(an2);
	mpf_clear(bn);
	mpf_clear(bn1);
	mpf_clear(bn2);
	mpf_clear(tn);
	mpf_clear(pn);
	mpf_clear(pi);

	// calcula o numero clocks necessario para executar o algoritmo
	t = clock() - t;

	printf("\nTempo de execucao paralela nova: %.3f s \n\n", (double) t / CLOCKS_PER_SEC);
	return 0;
}
